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The extrasolar planets found with radial velocity surveys have masses ranging 
from several Earth to several Jupiter masses. According to the conventional 
• sequential accretion scenario, these planets acquired super-earth cores prior to 

the onset of efficient gas accretion with rates that rapidly increase with their 
masses. In weak-line T-Tauri disks, mass accretion onto protoplanetary cores 
may eventually be quenched by a global depletion of gas, as in the case of Uranus 
and Neptune. However, such a mechanism is unlikely to have stalled the growth 
of some known planetary systems which contain relatively low-mass and close- 
in planets along with more massive and longer period companions. Here, we 
suggest a potential solution for this conundrum. In general, supersonic infall of 
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surrounding gas onto a protoplanet is only possible interior to both of its Bondi 
and Roche radii. At a critical mass, a protoplanet 's Bondi and Roche radii are 



X 



equal to the disk thickness. Above this mass, the protoplanets' tidal perturbation 
induces the formation of a gap. Although the disk gas may continue to diffuse 
into the gap, the azimuthal flux across the protoplanets' Roche lobe is quenched. 
Using two different schemes, we present the results of numerical simulations and 
analysis to show that the accretion rate increases rapidly with the ratio of the 
protoplanet's Roche to Bondi radii or equivalently to the disk thickness. In 
regions with low geometric aspect ratios, gas accretion is quenched with relatively 
low protoplanetary masses. This effect is important for determining the gas- 
giant planets' mass function, the distribution of their masses within multiple 
planet systems around solar type stars, and for suppressing the emergence of 
gas-giants around low mass stars. Finally, we find that the accretion rate into 
the protoplanets' Roche lobe declines gradually on a characteristic timescale of 
a few Myr. During this final stage, the protracted decline in the accretion flow 
across their Roche lobe and onto their protoplanetary disks may lead to the 
formation and retention of their regular satellites. 
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Introduction 



I n the convention al planet formation scenario ( ISafronov fc Zvjaginalll969l : iPollack et al. 



19961 : llda fc Linl 12004 ) heavy elements first condense into grains which then grow into plan- 
etesimals. In a gaseous environment, the hydrodynamic a nd tidal interaction between th e 
disk gas and these planetesimals damps their eccentricity (lAdachi et al.lll976l ; IWardlll998l ). 
Nevertheless, these planetesimals can perturb each ot her's orbit and induce growth through 
coagulation (lAarseth et al.l Il993l ; iPalmer et al.lll993l ). This growth process is stalled with 
the emergence of protoplanetary embryos that have acquired all the residual planetesimal s 
within their feeding zone, which is 3-5 times that of their Roche radius (IKokubo fc Idalll998l ). 
In protostellar disks with similar mass and temperature distributions as those inferred for 
the primordial solar nebula, the embryos atta i n isolation mas ses that are typically a few 
times that of the Moon at 1AU flLissauerlll987l : llda fc Linl 12004 ). 



In principle, the protoplanetary embryos can accrete gas from the disk when they have 
attained sufficient mass such that their surface escape speed is larger than the sound-speed 
in the disk. However, the gas which falls onto the embryos also conv erts gravitational energy 
into h eat which is then transported to the surface of the nascent disk (IBodenheimer &: Pollack 
19861 ). Unless the cores have acquired more than a few earth masses (M©), the envelopes 



aroun d them are too tenuous to allow efficient radiative or convective transport (IPollack et al. 



19961 ). Consequently, a quasi hydrostatic equilibrium is established and gas sediments onto 
the cores on a protracted cooling timescale. 

Once the protoplanet acquires sufficient mass to allow gas to be freely accreted, it enters 
a dynamical gas accretion phase. During this stage, the cores' gravity dominates the flow 
out to the Bondi radius (-Rb), where the local escape speed becomes comparable to the gas 
sound-speed in the disk. This radius increases linearly with the total mass of the protoplanet, 
which includes that of the core and the envelope. The protoplanets' Roche radius (Rr) also 
increases with the mass of the protoplanet, though the mass dependence is much weaker. 
Although Rb is generally smaller than Rr for low-mass protoplanets, Rb overtakes Rr for 
sufficiently large protoplanetary masses. 

As gas within R R collapses toward the cores, it is replenished by the diffusion of gas from 
beyond Rr. The rate of matter diffusion is generally assumed to be induced by turbulence, 
though the origin of such turbulence in protostellar disks at a few All's remains an out- 
standing issue. Protoplanets also exert a tidal perturbation on the disk with a strength that 
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increases with their mass. Around relatively massive protoplanets, this tidal perturbation is 
sufficiently strong to induce the for mation of a gap in the disk near the protoplanets' orbits 



( Goldreich &: Tremaind 1 19791 . Il980l ) and to quench the accretion (ILin fc Papaloizoul 11979 



19851 ). In relatively cold and weakly- viscous disks, the asymptotic mass of the protoplanets 



near the ice line in the solar nebula is comparable to that Jupiter (ILin fc Papaloizoul 11993 
Brvden et al1ll999h . 



However, there are numerical simulations of disk-planet interaction which in dicate that 
gas c o ntinues to flow in t o the gap after its form ation conditions have been satisfied (ILubow et al 
19991 ; iKley et al.ll200ll ; iD'Angelo et al.ll2002l ). In these simulations, the magnitude of both 
implicit and numerical viscosity is sufficiently large that the intrinsic redistribution of angular 
momentum in the disk gas is comparable to the tidal torque exerted by the embed ded proto- 
planet. Several simulations also t ake into account the effect of the turbulent flow (lArmitage 
19981 ; iNelson fc Papaloizoul 120031 ; lLaughlin et al.ll2004l ). All these models show that gas con- 
tinues to flow to the vicinity of the protoplanet until its asymptotic mass becomes several 
times that of the Jupiter. 

Although radial velocity surveys have discovered many planets with masses comparable 
to Saturn, the mass of Saturn is well below that found in the aforementioned simulations. 
While Saturn may have for med during th e late stages of protostellar evolution when the disk 
mass was severely reduced (lGuillotll2005l ). such a scenario cannot account for the relatively 
low mass of several close-in and resonant planets since it is generally accepted that their kine- 



mati c configuration requires extensive migration through active protostellar disks (ILin et al. 



19961 ). In this paper, we consider another effect which may limit the asymptotic mass of 
protoplanets involving the tidal barrier. Due to its multi-dimensional nature, the complex 
flow pattern is ideally analyzed with numerical simulations. 

The earliest disk simulations did not have sufficient resolution to resolve flow near the 



protop lanet (IMikil Il982l ). A novel approach was introduced by iKorycansky fc Papaloizou 
( 119961 ) to analyze inviscid flow in a steady-state. Under these conditions, the potential vor- 



ticity (hereafter vortensity) and the Bernoulli energy are conserved. In their simulations, they 
focused their attention on a local patch close to the protoplanet by adopting a set of boundary 
conditions which approximate the nearly co-orbiting flow at large azimuthal distance with a 
Keplerian velocity and smooth surface density. The advantage of this approach is that the 
high resolution allows the flow pattern to be accurately simulated (IBalmforth fc Korycansky 
200ll ). In the absence of viscosity, however, the stream lines around the protoplanets are not 
connected to that beyond the Roche radius. 



Adopting the same boundary condition, iTanigawa fc Watanabd (120021 ) carried out a 
series of high resolution time-dependent nonlinear simulations which clearly demonstrate the 
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presence of shock dissipation near the accreting protoplanet and generation/dissipation of 
vortensity across the shock. Nevertheless, between the shock waves, quasi steady flows can be 
established with conserved vortensity (tu) and Bernoulli energy Eb- They also demonstrate 
the possibility that a quasi steady-state may not be attainable under some circumstances. 
The minor deficiency of both these simulations was the artificial boundary condition at large 
azimuth from the protoplanet. We demonstrate here that the modification of the global disk 
structure can significantly alter the accretion rate onto the embedded protoplanets. 



There are many 2D global simulations of protoplanet-disk interaction (jPapaloizou fc Terquem 



20061 . see the recent review by). These simulations show that tidal torque from a sufficiently 
massive embedded protoplanet not only induces dissipative shock in the disk but also leads 



to gap form ation (jLin &: Papaloizoul 11993c lArtymowicz &: Lubowl Il994j ; iBryden et al.l [1999 



Kleyill999l . cf.). There are also 3D simulations that s how gas flow into the Roche lobe of 
the protoplanet (IBate et al.l 12003c iKlahr fc Kleyl 120061 ) . and some of these simulat ions have 
sufficient resolution to resolve the flow near the protoplanet ( D'Angelo et al.lbood ). The ac- 
cretion rate determined from these simulations suggests that despite a significant reduction 
in E near the orbit of the protoplanet, the growth timescale for Jupiter-mass protoplanets 
is shorter than the global-disk gas de pletion timescale (V ci m ~ 3 — 10 Myr) inferred from the 



obser vation of protoplanetary disks (IHaisch et al.l 12001 



Hartmann 2005; Silverstone et al. 



20061 ). If the gap region is continually replenished by diffusion into the gap, the asymptotic 
mass of t he protoplanets wo uld well exceed the observed upper limit of the planetary mass 
function jMarcv et al.ll2005h . 



Here, we explore the possibility that even if the gap is replenished with tenuous gas, 
only a small fraction of that gas may be accreted onto protoplanets with Roche radius (Rr) 
greater than the disk thickness H. The model parameters in many previous simulations are 
chosen to approximate those inferred from the minimum mass nebula model in which the 
ratio between the disk thickness and radius (H/a) is assumed to be around 0.07. In these 
simulations, H determines the magnitude of the sound-speed and therefore the pressure 
gradient in the radial direction in terms of the surface density distribution S. With this 
assumption, Rb becomes comparable to Rr for a protoplanet with M p ~ Mj. 

One exception to this set of parameters was the consideration of a pert urbation by a 
one-t hird Neptune-mass planet in very thin/cold disk with H/a ~ 0.02 — 0.03 (IBryden et al. 



20001 ) . In this case, Rb is several times Rr and the surface density of the disk diminishes well 
below that of an unperturbed disk. Although gas was not allowed to accreate onto the core 
of the protoplanet, the accretion timescales inferred from these simulations is substantially 
longer than those inferred from the idealized Bondi analysis, the viscous diffusion of the disk, 
and the Td ep inferred from the observed SED evolution of disks in young clusters. 
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The disk thickness (H) is determined by heating from both the intrinsic viscous dissi- 
pation and surface irradiation from the central star. The commonly assumed aspect ratio 
(H/a ~ 0.07) at the present-day distance of Jupiter (5.2AU) is determined under the as- 
sumption that the disk is heated to an equilibrium temperature by the proto-sun with the 
present-day luminosity of the Sun. Self consistent models indicate that the H I a at the snow 



line c an be significantly smaller than that assumed in the solar nebula model (IGaraud fc Lin 



20061 ). Taking this possibility into account, we suggest that, in the limit of Rr > H, a tidal 
barrier imposed by a protoplanet limits the gas flow from most regions inside the gap into 
the Roche lobe. Much of the mass that diffuses into the gap accumulates into the horseshoe 
orbits in the co-orbital region, but is not able to accrete onto the planet. 

In order to examine the dominant physical effects which regulate protoplanets' asymp- 
totic masses, we calculate the accretion rates and flow dynamics through two sets of sim- 
ulations. Since we are primarily in terested in the limit Rr > H, the flow pattern is well 



approximated by an 2D simulation (ID'Angelo et al.ll2002l ). In §2, we present an analysis to 



highlight the nature and importance of the tidal barrier. Due to the complexity of the 2D 
flow, we cannot obtain a full set of analytic solutions. For illustration purposes, we also 
adopt a grossly simplified ID approximation for the protoplanet's tidal potential and obtain 
solutions to highlight the effect of the tidal barrier in quenching the accretion flow onto it. 

In order to illustrate the structure of the flow, we present in §3, a set of high- resolution 
2D hydrodynamic simulations to demonstrate that the flow into the Roche lobe is indeed 
quenched for protoplanets with R R > H despite diffusion of gas into the gap region. Although 
the 2D simulation provides a more realistic approximation for the multi-dimensional geome- 
try of the potential field and the non-spherically symmetric flow pattern of the background 
gas, numerical simulations of this kind can only be carried out with a limited resolution for 
several dynamical timescales. However, the disk environment, as well as the mass growth and 
envelope structure of the protoplanet, all evolve on much longer timescales, especially during 
the epoch when the envelope undergoes a transition from quasi hydrostatic contraction to 
dynamical accretion. The computation of the late stages of giant planet formation also re- 
quires following the governing equations over large dynamical range in space and time. While 
the envelope around the core adjusts to quasi hydrostatic equilibrium, the flow around the 
planet's Hills radius adjusts to a steady state. The technical challenges for 2D simulations 
covering more than four orders of magnitude in physical scales over such large timescales 
makes them impractical for this purpose. In §3.31 we show that the dominant physical ef- 
fects generated by the 2D simulations can be quantitatively captured by a simple-to-use ID 
scheme, with which we verify that the gas accretion is significantly suppressed in the limit 
Rr > H. Utilizing this model, we show that the flow of gas into the protoplanets' Roche 
lobe declines and the protoplanets acquire asymptotic masses on Myr timescale. Finally, we 
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summarize our results and discuss their implications in §4. 



2. A Tidal Barrier for Accretion Flow 

In this section we begin our analysis with an analytic approximation of the accretion 
flow pattern in two- dimensions. Since protoplanets accreate from a disk environment, it is 
informative to compare the disk structure with the gravitational potential of the protoplanet. 
The scale height in the direction normal to the plane of the disk is determined by the gas 
sound-speed c s such that 

h = c s /n, (i) 

where Q = (GM^/a 3 ) 1 ^ 2 is the Keplerian frequency of the disk and M* is the mass of the 
central star. The Bondi radius for a steady, spherically symmetric accretion flow onto a 
protoplanet with a mass M p is 

R B = GM p /2c 2 s . (2) 
The Roche radius for a planet with semi major axis a is 



The magnitude of 



and, 



R R = (Mp/ZMjWa. (3) 



Rb ^Rr \3<?; 



(5) 



H _ 2h\ 
Rb Q3 

where q% = M P /10~ 3 M* and hi = c s /(0.1Qa) a re the normalized p lanetary mass and disk 
scale height. In a minimum mass nebula model ( Hayashi et al.|[l985l ) where the gas temper- 



ature is assumed to be identical to the equilibrium temperature of the large grains, hi ~ 0.7 
at the present radii of Jupiter (a = 5.2 AU). 

Under these conditions, the critical stage where the accretion rate is being suppressed 
occurs when R R ~ R B ~ H and M p = Mj. In principle, this region should be analyzed with 



full 3D consideration (IBate et al.l 120031 ; iKlahr fc Kleyl 120061 . cf.). However, the inadequacy 



of present-day facilities and the enormous computational requirements of full 3D numerical 
simulations place severe limits on both numerical resolution and temporal range. It is there- 
fore useful to first consider a more idealized 2D flow pattern. Through numerical simulations 
and analysis in this section, we emphasize that the protoplanets' accretion rate is sensitive 
to the global evolution of the disk structure. This approximation is also reasonable in the 
limit that Rr > H or equivalently Rr > Rr- 
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There have been multiple simulations studying the process of accretion onto a proto- 
planet covering a wide range of app roaches. In an attempt to achieve high resolution near the 
proto planet, several investigators (IKorycansky &: Papaloizoul Il996l ; iTanigawa &: Watanabe 



20021 ) considered flows in a local region surrounding the protoplanet. Both studies assume 



that gas entering the simulation from both the leading and trailing azimuthal boundaries 
has a constant surface density and Keplerian velocity. A second set of simulations have 
investigated the global disk response to protoplanet-disk interaction. Some of these simula- 
tions include the protoplanet's orbital migration while others include the effect of intrinsic 
turbulence. In all cases, they show that Jupiter mass protoplanets induce gap formation in 
protostellar disks with properties similar to that of the minimum mass nebula. In order to 
achieve sufficient resolution for studying the flow around protoplanets, most of these simu- 
lations are carried out for Mj-mass protoplanets in disks with H/a ~ 0.07 such that their 
Rr ~ Rb ~ H. Here, we are primarily interested in the disk regions where Rb » Rr 
which is possible only for Rr > H. We show that gas-giants with M p < Mj/3 (such as 
Saturn) are probably formed in regions of the disk that satisfy H/a < 0.04 during the epoch 
of planet formation. 



2.1. Streamlines Intersecting The Feeding Zone 

The continuity and momentum equations for a two-dimensional inviscid disk can be 
obtained through vertical integration such that 

— + V-(Sv) = (6) 



— + v • Vv = -2fizxv - V$ 9 - c^VlnS. (7) 



and, 

~dt 

Similar to the ID analysis presented in §2.2[ we have adopted an isothermal equation of state. 
We determine the magnitude of v in a frame which rotates at the same angular frequency Q 
as the planet orbits its host star. In the analytic derivations, we neglect the effect of viscous 
and shock dissipation, but some degree of dissipation is present in our numerical simulations. 
These approximations are adequate for linking streamlines between shock waves. We also 
adopt a Hill's approximation for the potential in this frame such that 

$ =--n 2 x 2 ^ (8) 



where the shearing- sheet coordinates (IGoldreich fc Lynden-Belll 119651 ) are centered around 



the planet with x and y correspond to the radial and azimuthal directions. Although this 
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expression is a much more accurate description of $ 9 than the one- dimensional approximation 
given in Equation (fT9]) . approximations must be made to solve it analytically. 



In the absence of any viscous or nonlinear dissipation, it is helpful to analyze the flow 
with a vortical stream function (ip) and a scalar potential (4>) which varies along the lines of 
constant ij) such that the 2D mass flux can be decomposed into 



Ev = Vxz^ + V0. 



(9) 



Contours lines of constant ip and are orthogonal to each other, and a steady-state (V ■ 
(Sv) = 0) corresponds to V 2 = 0. The velocity is parallel to streamlines of constant x/j, 
along which the mass flux (Su) is conserved, and the momentum equation can be reduced 
to two integrals of motion 

^ = do) 



v 2 /2 + c>£ + $ 9 = E B (1>) 



(11) 



which correspond to the vortensity ( 


PaDaloizou & Lin 


1989 


) and Bernoulli constant ( 


Korvcanskv & Papaloi 


1996 




Balmforth & Korvcanskv 


2001 


)■ 







Differentiating the mass flux F m = Ev and Eb with respect to <fi along these streamlines, 



(dhi£/<i0)^ = — (dliav^/dcf))^, 



and, 



\d<t>). 



dx 



dy 



x 



3Q 2 



GM„ 



(x 2 + y 2 fl 2 



GM„ 



T,Vy J (x 2 + y 2 ) 3 / 2 



(12) 

(13) 
(14) 



where we define the magnitude of v along the streamline to be ity. Equation ( fl4T) clearly 
indicate a potential transonic point along each streamline. 



2.2. Analogous Bondi Solution 

The sonic transition in a 2D flow is similar to that in the ID flow, and the quantities 
F m and Eb are equivalent to the mass flux 

M p = Airpur 2 (15) 
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and the one-dimensional Bernoulli energy 

E m = u 2 /2 + c>p + $ 9 (16) 

where p and w are the spherically symmetric density and radial velocity. We find it is useful 
to identify the dominant physical effects which may affect the 2D flow pattern with the help 
of an approximate ID solution. 

In a steady state, the equations of ID continuity and motion are combined into 

(c 2 -u 2 )du 2c 2 „ x .„„. 

IJ L— = 2--V$. (17) 

u dr r 

In the standard Bondi analysis for a point-mass potential in which $ 5 = $ = —GM p /r, the 
continuity of a steady flow across the sonic point r s where u = c s requires r s = —2c 2 /VQ 
and Eb = c^lnpoo. By scaling the background quantities with those at the sonic point r s 
such that p s = pooexp3/2 ~ 4.5poo, we find the accretion rate to be 

M p = M pm ~ U Poo G 2 M 2 /c 3 s . (18) 

In this case, the potential at r s is lower (more negative) than at r = oo, and the gravitational 
focusing effect causes p s > p^. 

However, around a protoplanet with a circular orbit and semi major axis a, the region 
outside its Roche radius (Rr = (M p /3M*y/ 3 a) is dominated by the gravity of its host star. 
In the frame which co-rotates with its orbit, the full tidal potential due to both the planet 
and its host star is given by Equation (jSJ). For the ID analysis in this section, we introduce an 
approximation to the gravitational potential which includes the tidal effect induced by central 
star but neglects its gravitational torque on the otherwise assumed spherically symmetric 
(relative to the protoplanet) radial flow such that 

*^*°H# + itr)- « i9) 

Here we did not include the Coriolis force associated with the rotation of the frame. In this 
approximation, gravity reduces to that of a point mass (V$ — GM p /r 2 ) for r < Rr, it be- 
comes repulsive, (V$ ~ -GM P /[R 2 R (1 + r / R^j + GMp/r 2 ) for R R < r < R^, and vanishes 
asymptotically (V$ — —GMpR^/R^r 2 ) for r >> R^. In principle, this approximation can- 
not be used to replace the full mult i- dimensional simulation on the detailed structure and 
dynamics of the outer most region of the envelope near the protoplanet 's Roche radius. How- 



ever, by setting R^ = r max = \/12Rr (IGreenzweig k, Lissauerl Il990l ). the potential in the 



frame which co-rotates with the protoplanet's orbit can be reasonably approximated along 
the radial direction. 
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For the tidal potential, the right hand side of Equation ffTTj) vanishes at r s = Rr/£ 
where the parameter £ can be obtained from a quadratic equation 

^ \ Roc Rb) ^ \R-lo RodRb ) RbRLC, ^ ^ 

the solution of which is plotted in Figure In the limit that Rr < Rb, £ — 1 — Rr/Roo + 
Rr/2R b , i.e. the sonic point is close to the Roche radius. 
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Fig. 1. — Sonic radius parameter £, defined by the sonic point r s = Rb/C,, as a function of 
the ratio of Roche to classical Bondi radii. The dependence of £ on R R /R B can be found in 
Equation (J2D). 
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When the tidal potential is taken into account, Er = c^lnpoo — GMpR^/ R? R at large r, 



such that the density at the sonic point is given by 

~$(r s ) - $fool 

Ps = Pooexp 



(21) 



Since the absolute value of the potential at r s is near a maximum (where gravity V$ ~ 0), 
Ps < Poo- This density decline is due to the need of a positive pressure gradient, required to 
overcome the tidal barrier associated with this potential maximum. In the limit of relatively 
small c s (or equivalently large Rb) a large density reduction is needed. Based on the above 
extrapolation that r s ~ Rr for Rr << Rb, 



where rj = 2(R OD /R R - 2) 
r s ~ R R , and 

2tt 



Ps ~ 0.6 Poo exp (-tjRb/Rr) (22) 
3 is positive. In a tidal potential with Rr « Rb, £ ~ 1, 



M p = M t ~ ( c^pooexp (^7^ 



0.135M 



exp l — 



(23) 



The above equation indicates that, in the limit Rr « Rb, the tidal accretion rate (M t ) can 
be much smaller than the standard Bondi accretion rate for a point mass potential M pm . 



2.3. Boundary Conditions in 2D Flows: Loading the Streamlines 

A comparison of the governing equations (1171) and (TH|) indicates that the competing 
effects of pressure gradient and gravity in the 2D flows is determined by the boundary 
conditions of the streamlines. At large x distance from the planet, the unperturbed 2D 
flow is approximately Keplerian, with v x ~ and v y ~ — 3Qx/2, which can be substantially 
larger than the sound-speed. These streamlines are on hyperbolic paths with large impact 
parameters relative to the protoplanet and material on these streamlines will not accreate 
onto the planet. 

The streamlines which are relevant to accretion flow are those with impact parameters 
smaller than Rr. In the limit H « Rr, pressure does not contribute significantly, and 
Equation ([7]) reduces to that for non-interacting particles. Particles that originate from un- 
perturbed Keplerian orbits at = x max = \^12Rr and i/oo = ±ica, reach the vicinity of the 



unsta ble Lagrangian saddle points (±Rr, 0) with negligible velocity (IGreenzweig fc Lissauer 



19901 ). The region of the disk between a — x^ and a + x^ is commonly referred to as the 



feeding zone. Although particles with initial Keplerian orbits that pass through the location 
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Uoo = ±7ra with Xoo = 1.8a; — 2.7 Rr can directly strike a co mpact protoplanet co re, particles 



on neighboring orbits are strongly scattered by the planet (IBryden et al.ll2000l ) 



Although w and Eb are useful quantities for analyzing the flow in the horseshoe region, 
they are not conserved for gas accreted from the disk onto the planet. More specifically, for 
fluid in the protostellar disk with Keplerian asymptotic velocities, w = fi/2£oo. However, if 
gas enters into the planet's Roche radius with x 2 +y 2 < 3R 2 R /A in a Keplerian circumplanetary 
disk, its streamlines have uo < —2Q and w < 0. This mismatch in w implies that gas can 
only be accreted from protost ellar disks onto protoplanetary disks through shock dissipation 
(ITanigawa fc Watanabdl2002l . see 2D simulations by). Since any dissipation would invalidate 
the conservation of w and Eb, the flow pattern derived from th ese conservative quantities 
cann ot adequately describe accretion flow onto the protoplanet (IKorycansky fc Papaloizou 
19961. cf.V 



However, these quantities are useful for determining the conditions under which gas 
on some streamlines in the horseshoe regions may avoid shock dissipation and accrete onto 
the protoplanet. The flow in this region depends sensitively on the boundary condition at 
t/oo). If the flow at (xoo,2/oo) is strictly Keplerian, as assumed in some previous local 



simulations, Eb = cjlnEoo — 3^2 x^/8. In this case 



S s = Yj(x s ,y s ) — Sooexp 
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H 2 2 



(24) 



at the sonic point (x s ,y s ) where u = c s . The parameter 3 is independent of H or c s . For 
Keplerian boundary conditions, the value of 



xz 



R R 



2R r 



X 



(25) 



In this expression for 3, the actual values of (x s , y s ) can only be determined by the detailed 
solutions of Equations (jSJ) and ([7j). 

In analogy with the ID Bondi solutions, we suggest that there may exist a sonic tran- 
sition point in Equation (fl4l) as is found in Equation (ITT]) . In this case, if \u\ < c s at some 
point along the streamline, the right hand side of Equation (1T41) would vanish during the 
sonic transition at (x s ,y s ) with 



Vx(x s ,y s ) = Xs( (x 2 s + Vs) 3/2 



(26) 



Vy{x s ,y s ) 

The point (x s ,y a ) is close to the turn-around point (where u y vanishes) of the horseshoe 
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orbit for non-interacting particles. Since v 2 + v 2 — c 2 at (x s ,y s ), we can express, 

-i 



(27) 



In the ID approximations (see the previous subsection), the physical parameters at 
the sonic points can be analytically determined from those at infinity But in 2D flows, 
Equations ( fT4l) and ( !27l) generally do not lead to an analytic relation between r s and c s . 
However, their approximate solutions suggest that the condition for v 2 < c 2 and v 2 < c 2 is 
r 2 = x 2 + y 2 ~ R 2 R) i.e. in the 2D limit, the transonic radius is also comparable the planet's 
Roche radius (r s ~ Rr) and 



M p ~ 2irR R c s E s = 2TcR R HttY loo exp 



\H 2 



(28) 



The above expression clearly indicates the importance of (3. I f /3 < 0, M will decrease with 
H as shown in the local simulations of Tanigawa fc Watanabel (2002). However , if (3 > 0, M 
will increase with H as shown in the global simulation by iBryden et al.l (119991 ) . 



We suggest this dichotomy can be attributed to the boundary conditions at (x^, y^). In 
order to illustrate this conjecture, we first consider a special streamline which passes through 
a sonic transition at the same azimuthal phase as the planet (i.e. y = 0), the regularity 
condition for Equation (fl~4l) requires d§ g /dx = d§ g /dy = or equivalently, x s = Rr. In the 



limit that H « R R , $ g (x s ,y s ) = -9Q 2 R 2 R /2, x c 



Vl2R 



R, Vo 



ira, and (3 ~ 0. 



These streamlines are analogous to the parabolic trajectories of those particles which can 
marginally reach the separatrix at the Lagrangian points from an initial circular Keplerian 
orbit. Starting with Keplerian motion, streamlines with > x max do not enter into the 
planet's Roche lobe. However, Keplerian streamlines with R R < x^ < s max can enter into 
the planet's Rr and cross the sonic radius with very small negative values of j3. In local 
simulations done previously, the disk flow at some specified boundaries is explicitly set to 
be Keplerian at an azimuthal phase y^ < ixa. This artificial choice of boundary condition 
leads a small negative value of (3 for disks with both negligible and finite H J s. This result 
implies that tidal barrier is overcome purely by the initial kinetic energy rather than the 
positive pressure gradient and T, s is slightly larger than as a consequence of convergent 
flow. If (3 is negative, Equation (|2"5I) in dicates that M p is a dec r easing function of H. This 
is in agreement with the simulations of iTanigawa fc Watanabd (120021 ) . With their value of 
Soo (~ S mn ), the growth timescale would be r p = M p /M p ~ 10 3 yr(M p /10M e )" a3 , which 
is much shorter then the global depletion timescale (r^ ep ) for the disk. The discrepancy 
between these timescales illustrates the need for a mechanism to quench the accretion flow. 
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In contrast to local simulations, full-scale 2D simulations the protoplanet induces the 
formatio n of a gap in the disk through both accretion onto itself and tidal dissipation within 



the disk (IBryden et al.lll999l ). In such simulations the boundary conditions at (x OQ ,y 



are 



significantly different from those assumed in local simulations. Near the edge of the gap, 
the perturbed pressure gradients cause the velocity of the gas to become super Keplerian at 
positive Zoo and sub Keplerian at negative Xoo. The magnitude of the flow velocity, 



Vy \Xc 



Hoc) = Via 



H 2 <91n£ 
2a 2 <91nx 



3x\ 



(29) 



is reduced for both positive and negative value of x^. The typical gap has a width R R < 
A < x max . When the thermal and viscous conditions for gap formation (Rr > H and 
Mp/M* > 40v/tta 2 where v is the effective viscosity) a re satisfied, the surf ace density inside 
the gap £00(^00,2/00) ~ eS with e < 1(T 2 - 10" 4 farvden et al.l Il999h . Although this 



reduction in decreases the value of M p (by a factor similar to e), the growth timescale 
t p would still be more than an order of magnitude shorter than Td ep if the protoplanet 
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20031 ) or due to intrinsic or artificial viscosity (lArtymowicz fc Lubow 



In the limit of severe surface density depletion near the protoplanet, the magnitude of 
within the gap may be substantially reduced from the Keplerian approximation. In the 
limit \xno\ < H, we approximate Equation ([29]) as v?,(xoa, Uoo) = (9/4)fi 2 x 2 XD (l — / su b) where 

/sub = (30) 



< H, we approximate Equation (1291) as v^x^ 

2HH dln£ 
3 x a dlnx 



For values of < / < 1, the Bernoulli constant can be written as Eb — cpnEoo — 3(1 + 
3/sub)^ 2 ^L/8- At large azimuthal phases, the streamlines (constant ip) are nearly parallel to 
the y-axis. Along these streamlines, the potential has a maximum at d& g /dx = 0, which is 
located on a circle x 2 + y 2 = R 2 R centered on the protoplanet. In order for the background gas 
to reach the protoplanet within the Roche radius, it must pass through this local potential 
maximum with a positive pressure gradient. In the limit of low sound-speed, a sufficiently 
strong pressure gradient is only attainable with a large density gradient, i.e. the severe gas 
depletion interior to Rr. For these conditions, the value of 



(3 = 3 



1 {( 1 + 3^x1-4x1 

R R 



(31) 



In the limit x c 



Xr\ 



and x s ~ Rr, (3 ~ 27/ su b/2 > 0. For any value of < / su b < 1, (3 > 



for all streamlines with x^ > (vT2/(l + 3/ su b)) Rr. Streamlines with smaller values of x^ 
attain horseshoe orbits that do not enter into the Roche radius. In cold disks (H < Rr) this 
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tidal barrier suppresses the flow in the azimuthal direction such that only a small fraction of 
the gas which is diffused into the gap can actually reach its Roche lobe. Consequently, the 
accretion rate derived from Equation (I28p become substantially smaller than that obtained 
when neglecting the effect of the azimuthal tidal barrier. 



2.4. Protoplanet's Asymptotic Mass 

We now apply the appropriate boundary condition in protostellar disks to determine 
the asymptotic mass of protoplanets. For illustration simplicity, we first adopt the ID 
approximation and then consider the full 2D flow. 

The density of the ambient gas can be derived from = T,/2H where the surface 
density is £ = f g T, mn and f g is a scaling factor. We use the minimum mass nebula model 
E mn ~ E (a/1AU)- 3 / 2 , where S (t ) = 2 x 10 3 g/cm 3 fid ucial prescriptio n. We also 



approximate the evolution of the disk by E = S exp(— t/r^) (jlda fc Linll2004j ) where the 
characteristic gas depletion timescale Td ep inferred from the observed SED is t& cv ~ 3 — 
WMyr. 

From the ID mass accretion rate given in Equation f[2"3"j) . we can determine the asymp- 
totic mass of the planets Mf by first solving for Xf — 3 1 ^(r]/2) 1 ^ 2 h^ 1 q^j 3 with 

Xf expx'dx = A l^^^ (32) 

where A = (10 8 2 3 /3)^ 2 tc 2 !] 1 ^ 2 is a constant, and the planet's asymptotic mass is contained 
in q 3 f = Mf /10~ 3 M*. The left hand side of Equation ( 132|) is a rapid increasing function of Xf. 
For a minimum mass nebula model [i.e. f g — 1 at a=5 AU) the magnitude of its right-hand 
side is ~ 10 6 and the solution of Equation (1321 is Xf ~ 4. However, the protoplanet may be 
embedded inside a gap where the gas density is reduced by at least 2 orders of magnitude 
from that in the minimum mass nebula. Long period protoplanetary embryos may also have 
formed during the advanced stages of disk evolution when the magnitude of f g may also be 
relatively small. However, the value for Xf is relatively insensitive to the magnitude of f g 
and it reduces to ~ 3 if f g ~ 10 -2 . The corresponding value for the asymptotic planet mass 
is 



The structural parameters of the disk are functions of M* and a. Near the snow line of a 
minimum mass nebula, the magnitude of Mf is comparable to that of Jupiter. Equation fl33l) 
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is also equivalent to Rn(Mf) ~ H. This relation is in agreement with the thermal condition 
for the tidal truncation and gap formation of protostellar disks by embedded protoplanets 
(ILin &: Papaloizoulll986l ). 



The dependence of S s /Eoo on the sign of (3 in 2D is similar to that of p s /poo on the 
presence of the tidal barrier in the ID approximation. Using the 2D mass accretion rate given 
by Equation (|28|) . together with the approximate global £ depletion formula, we deduce an 
asymptotic mass for a protoplanet of 



3 (H_ 



ln- 



57T 



HM*P 



3/2 



(34) 



where M gap = 7reE a 2 is the characteristic mass associate with the gap region. For /3 ~ 1 and 
e ~ 10~ 3 , the above asymptotic mass would yield Rr ~ iJ. However, if the magnitude of j3 is 
very small or negative, the above equation reduces to M p /M* ~ 3(87rifM gap rdcp/97raM*P) 3 / 2 , 
which is substantially larger. To calculate the actual values of j3 and f su b, we must utilize 
numerical simulations 



3. Two Dimensional Numerical Simulations 

3.1. Numerical Method 

In order to verify the 2D analytic approximation in the previous section and deduce a 
value for (3, we now present a series of 2D numerical simulations. The 2D numerical scheme 
we use is a fully parallel hydrodynamical code with which the continuity and momentum 
equations are solved on a fixed Eulerian grid in 2D cylindrical coordinates (r,<f>). A stag- 
gered grid is introduced where scalars are defined at grid centers and vectors are defined on 
grid edges, allowing us to maintain second-order spatial accuracy. Earlier versions of the 
code were used, among o ther applications, to study the formation of gaps in protostellar 



disks (IBryden et al . 2000) and also boundary layer physics in the accretion disks around 



cataclysmic binaries (lKleylll989l ). The equations of continuity and motion for the fluid and 
the expression for the gravitational potential are similar to those given in Equations (jBl), ©, 
and (jU) respectively. In the equation of motion, we adopt an isothermal equation of state 
P = c^p, with a radially varying sound-speed, c s = (H/a)Vk ep where the aspect ratio of the 
disk H/a is assumed to remain constant. In order to make a direct comparison with the re- 
sults of the analytic approximation, we neglect the effect of the viscous stress associated with 
the turbulent motion. A softening parameter r so f t = 0.5Rr is adopted for the gravitational 
potential of the planet. 
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The planet is initialized with a lMj up planet at a location (r, 0) = (1AU, 0). Although 
the planet is held fixed, any mass accreted during the simulation adds to the overall gravita- 
tional potential of the planet. The simulations cover the entire azimuthal extent of the disk 
and a radial range from 0.5AU to 1.5AU. Simulations were initialized with a surface density 



profile of £ = 950^3- (jjfj) 32 , and have a numerical resolution of 



(200, 1600). 



Each of the terms in the fluid equations are dealt with consecutively through both 
operator and directional splitting. Our advection scheme is an e xtension of the simple first- 
order upwind scheme based on the prescription of Ivan Leer! (119771 ). and allow s for the detailed 
treatm ent of strong shocks. Details of the advection scheme can be found in iKley fc Hensler 
(119871 ). We solve the entire problem in a frame rotating with the planet. 



3.2. Two Dimensional Numerical Models 



Before the fluid is incorporated onto the planet it must navigate though both the shocks 
and sub-disk surrounding the protoplanet. In the 2D study we neglect any radiative feedback 
onto the gaseous disk from the e nergy released during th e cool ing; and contraction of the 
planetary atmosphere. Following ITanigawa fc Watanabei (120021 ). we decrease the surface 
density of the disk in the region r sink < 0.5Rr at a rate given by 



E (t + At) = E (t) 



At 

Tsink 



(35) 



This mass is a dded to the mass of the protop lanet, and gradually increases its gravitational 
contribution. ITanigawa &: Watanabei (120021 ) demonstrated that the overall accretion rate 
(M) was relatively unaffected by the choice of r sink provided that r sink « 0.1R R . Because of 
our lower numerical resolution near the protoplanet we set r sink = 0.5Rr and r S i n k = 1/^1. 
The dominate effect which regulates the accretion rate is the supply of material to the sub- 
disk region. This quantity is determined by the effectiveness of loading the streamlines far 
from the planet, so the values of r sink and r S i nk should not be crucial in determining the 
asymptotic mass of the planets. In order to test the validity of this assumption and the 
convergence of our results, we carried out several simulations with varying values of T s i n k, 
but these showed little variation in the overall accretion rate. 

We present the results of two models in which H/a = 0.04 and 0.07 respectively. Figure 
(J2J) shows both the total mass and accretion timescale as a function of orbital period for these 
models. The accretion rates approach a steady state on the time scale of 100 orbital periods. 
This time scale is comparable to the liberation period of flows in the horseshoe region and 
the synodic period of the flow near the separatrix at the Lagrangian points. It reflects the 
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protracted adjustment of the disk structure to the tidal perturbation of the protoplanet. 
The accretion timescale (tiau = M p / M p )iaxj f° r a Jupiter- mass protoplanet is proportional 
to both the local surface density (oc a 3 / 2 ) and to its orbital period (oc a 3//2 ). It can therefore 
be scaled to other radii such that 

Tgrowth = riAu(a/lAU) 3 . (36) 

Figure ([2]) shows timescales for planet located at both 1 and 5AU. For a planet at 5AU, 
we see that the growth timescale becomes comparable to the gas depletion time scale of the 
disk. Within the computational duration, the total fractional gain in the planetary mass 
is 0.05 and 0.2 for H/a = 0.04 and 0.07 at 1 AU and one order of magnitude smaller at 5 
AU. Protracted 2D simulations for protoplanets' nonlinear growth (over 10 4 orbits) require 
computational resources beyond that currently available. 

The results in Figure (j2J) also clearly indicate that the mass of a planet embedded in a 
disk with a lower aspect ratio levels off at significantly lower mass than one embedded in a 
disk with larger (H/a). Based on these results, we compute, from Equation (1281) . the value 

where the subscripts 1, 2 refer to the values for the two models. The value of f3 increases 
from a very small value to the order of unity after ~ 50 orbits (see Figure (jBJ)). The disk is 
initialized as a Keplerian disk, and as a result there is some amount of mass on Keplerian 
orbits within the separatrix region (i?R < x < \/Y2Rrj at azimuthal phase < na. This 
gas passes over the tidal barrier near the planet solely based on this artificial kinetic energy. 
For this reason the first ~ 30 orbits in Figure ([2]) show similar accretion rates onto the planets 
in both simulations and the initial value of (3 is more closely approximated by Equation (I25p . 
In fact, during this relaxation phase, the lower sound-speed simulation is able to accrete mass 
onto the planet faster because the outward pressure gradient that develops around the planet 
is smaller. 

After approximately one libration period around the L 4 and L 5 points, mass elements 
that start their trajectory at the superior conjunction (i.e. y = an) can reach the vicinity of 
the protoplanet 's Roche radius. Thereafter, a large pressure gradient along the stream line 
is needed to overcome the tidal potential barrier. During this second phase, the accretion 
timescale can be up to an order of magnitude larger in the high sound-speed case. The value 
of (5 ~ 1 is more closely approximated by Equation (1HT1) . The discussion near Equation (1291) 
indicates that this approximation is based on the assumption that v y (x oc , y^) ~ 0. In Figure 
(jl]), we plot the value of v y at the planet's superior conjunction as a function of the disk 
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radius. We note that within the feeding zone, the magnitude of v y is significantly reduced 
from its Keplerian values which is in agreement with Equation (1291) . 
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Fig. 2. — The total mass of the planet as a function of time (top). The solid line shows 
results for (^) = 0.04, while the dotted line shows (^) = 0.07. The lower panel shows the 
accretion timescale for both simulations. The left-hand ordinate shows the timescale for a 
planet at 1 AU, while the right-hand ordinate shows the timescale for a planet located at 5 
AU. 
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Fig. 3. — The evolution of /3, deduced from the comparison of two models with f —J = 0.04 
and 0.07. From Equation ( |28i) . we see that a positive f3- value indicates that M is an increasing 
function of the disk scale height (if). 
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Fig. 4. — The radial distribution of v y at the protoplanet's superior conjunction for models 
with (^) = 0.04 (solid line) and 0.07 (dashed line). A decline in the magnitude of v y within 
the feeding zone is due to a pressure gradient associated with the surface density distribution 
across the gap. 
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Our simulations do not include explicit viscosity, and as a result mass diffusion into 
the gap is due to numerical viscosity or shocks. Figure (j3J) shows the total mass within the 
separatrix region (\r — a\ < ^/12Rr) for the two simulations. Again, during the initial stages 
of the simulation the two cases track each other. The drop in gap mass is quite precipitous 
(~ 0.3M Jup ) during the first 100 orbits. Comparison with Figure © shows that a majority 
of this gas is not accreted by the planet, but rather shocked by the planet and pushed out 
of the region. However, the lack of material is not the cause for the mass of the protoplanet 
to level off in the low sound-speed model; by the end of the simulation, there is more mass 
within the gap region in this simulation than in the high sound-speed simulation. 
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Fig. 5. — The total mass of the planet in the gap region (|r — a\ < \/Y2Rr} as a function 
of orbits for — = 0.04 (solid line) and — = 0.07 (dashed line). By the end of the simulation 
there is more mass within the gap for the low sound-speed case, indicating that the tidal 
barrier, not the lack of material within the feeding zone, is the cause of the planets increasing 
growth timescale. 
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In order to further illustrate that the pressure gradient is insufficient to overcome the 
tidal potential barrier along the stream line in the low sound-speed simulation, we identify 
streamlines with contours of constant vortensity w. Since shocks occur only relatively close 
to the protoplanet, vortensity is approximately conserved along the stream lines in regions 
outside the Roche lobe. In Figures © and (JTj) we overlay the contours of w on the velocity 
field. Outside Rr, the streamlines are nearly parallel, and those that eventually enter the 
protoplanet's Roche lobe are initially located at ~ 2 — 3Rr ~ 0.15 — 0.2a from its orbit. 

Having determined that the flow is primarily along contours of constant w, we can 
utilize the discussion at the end of §2.3 to link the flow velocity at the superior conjunction 
to the resulting accretion rate. In Figures (jHJ) and 0, we plot A = {y y + Qa) /v kep — 1, 
where v kep = (GM^/a) 1 ^ 2 and v y is given in the rotating frame. Comparing to Equation 
( 1301) . we see that / su b — 2 |A|. At a = 0.8AU and 1.25AU (which correspond to ~ ±x max 
away from the protoplanet), |A| ~ 0.03 and thus / su b ~ 0.06. For a streamline that connects 
%oo = imax and x s = Rr, Equation (|3T!) gives (3 ~ 0.8 — 1.2, which is in agreement with the 
results shown in Figure 
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Fig. 6. — Potential vorticity (w) contours overlying the velocity field of the disk for the low 
sound-speed simulation with (— ) = 0.04. Potential vorticity, conserved in regions outside 
Rr, provides a good tracer of the streamlines within the flow. 
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Fig. 7. — The same as Figure (EJ), but for (^) = 0.07. 
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Fig. 8. — A = (v y + f2a)/t> kcp — 1 at the protoplanet's superior conjunction for (^) = 0.04. 
This quantity parameterizes the deviation of the flow from Keplerian (A = 0). A is related 
to / su b — 2 |A|, and can be used to determine the magnitude and sign of (3, the accretion 
rate parameter. 
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In Figures ffTUl) and ffTTj) . we plot the variations of E along the azimuth at constant 
semimajor axis. Although the approximation of Equation (|14p breaks down near the pro- 
toplanet, it provides a reasonable indication of the variation of E along the streamlines. 
The results in Figures ( ITU]) and (ITTj) show much larger azimuthal amplitude variations in E 
for the low sound-speed model than for the high sound-speed model. The magnitude of E 
decreases by two orders of magnitude in H/a = 0.04 model. This result is consistent with 
our analytic discussion that a large density gradient is needed in the low sound-speed limit 
in order to sustain sufficient pressure gradient to overcome the tidal potential barrier, and 
can be compared to the ID results of Figure (fT5j) . 
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Fig. 10. — Surface density (£) as a function of azimuthal angle (0) at several values of 
constant r around the planet for f— J = 0.04 at i = 500 orbits (left panel). Contours are 
at r = 0.95, 1.0, and 1.05AU, while the planet is located at 1AU. The precipitous drop in 
density along azimuth provides the necessary pressure gradient force needed to overcome the 
tidal barrier. The right-hand panel shows the overall surface density in the region. 
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Fig. 11. — Same as Figure ( ITUl) but for f~) = 0.07. Given the higher sound speed, the 
density drop needed to overcome the tidal barrier is not as large as in the low sound-speed 
simulation. 
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Figures ffl~2]) and ffl~3l) show the Mach number ((^+^) 1,/2 /c s ) distribution in the rotating 
frame near the planet for (H/a) = 0.04 and 0.07 respectively. Both the approximately 
Keplerian flow far from the planet, and t he horseshoe orbits within th e feeding zone are 
clearly visible in both plots. As observed in lTanigawa fc Watanabd (120021 ). the radial extent 
of the shocks from the planet is proportional to the soundspeed. For the most part, the 
subsonic region is bounded by a strong shock. However, for the high sound speed case, 
the subsonic region also extends from inside the Roche lobe to the feeding zone beyond t\r, 
allowing gas to efficiently flow into the vicinity of the protoplanet. For the low sound-speed 
model, The subsonic region near the protoplanet is well separated from the low- velocity horse- 
shoe region. These figures once again support the conjecture that large density gradient is 
needed for the pressure to overcome the tidal potential barrier. 
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Fig. 12.— The Mach number for (f ) = 0.04 in a frame rotating with the planet. The 
dashed circle over-plotted on the flow denotes the approximate Hill sphere. The subsonic 
region surrounding the planet is well separated from the rest of the flow, and a large density 
gradient is needed to overcome the tidal barrier. 
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Fig. 13. — The same as Figure (pi) but for (f ) = 0.07. The subsonic region surrounding 
the planet is connected, via subsonic streamlines, to the material in the feeding zone. Thus, 
the density gradient needed to overcome the tidal barrier is not as large as in the (f ) = 0.04 
simulation. 
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3.3. Simulations Based on the ID Tidal Barrier Approximation 

The main objective in our simulation is to determine the asymptotic mass of protoplanets 
in depleting protostellar disks. The results of our analysis and numerical simulations (Figure 
(J2J)) indicates that protoplanets reach their asymptotic mass when their growth timescale 
becomes comparable to that of the disk depletion timescale, which is 10 5 time longer than 
their orbital period. The accreted gas collapses to form the envelope of gas-giant planets with 
size scales several times the present-day radius of Jupiter. The dynamical timescale near the 
core is ~ 10~ 5 that of the orbital period. While we would like to simultaneously study the 
evolution of the gaseous envelope and the decline of the mass influx through the protoplanet's 
Roche lobe, the limitation of currently available computational facilities prohibits the high- 
resolution, multi-dimensional simulations of the planetary accretion process over such a large 
range in the relevant timescales. 

Here we carry out a series of ID simulations based on the spherically symmetric ap- 
proximation of the tidal potential ( §2.21) . The basic equations for the ID code include three 
conservation laws of momentum, mass and energy (including radiative diffusion and con- 
vection), and the equation of state. During the phase of dynamical accretion , conventional 



numerical schemes based on a hydrostatic approximation (jPollack et al.lll996l ) are no longer 
adequate for the determination of the planetary structure and evolution. We developed a nu- 
merical scheme which is based on a Lagrangian approach where the radial coordinate evolves 
in time. It was developed to study the envelope structure as it undergoes a transi tion from 



quasi hydrostatic equilibrium to dynamical collapse (see a detailed description in (ILi &: Lin 



20061 )). The advantage of this method is that we can use it to resolve the large spatial 
range of the protoplanet's envelope and to follow its evolution over an extended duration of 
time. The main limitations of this method are the assumption of spherical symmetry and 
our neglection of the protoplanet's tidal torque on the global structure of the disk gas. We 
utilize the versatility of this scheme to highlight the dominant physical effects and compare 
our results to the 2D scheme to assess the validity of our ID calculation at various stages of 
protoplanetary growth. 

For the present application, gas is accreted onto a Saturn-mass object (M p = 100M e ). 
During the dynamic gas accretion stage of planetary formation, the envelope of the proto- 
planet will contract much more rapidly than before. While the gas flow close to the core 
remains in quasi-equilibrium, the outer region of the envelope, especially the region near 
the surface of the proto-planet, should experience dynamic gas accretion. We focused our 
calculations on the dynamic region of the gas flow, which covers the radial region from 10 core 
radii to 2 Bondi radii. Across the inner boundary, we calculated the gas accretion rate and 
infer a luminosity based on the assumption that the inflow is halted at the core radius. As 
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the gas flow in this region tends to a steady-state (to be shown below), the radiative transfer 
is very efficient and there is insufficient thermal pressure to counter balance the core's gravity 
throughout this region. It is therefore a reasonable assumption for the velocity across the 
inner boundary to be close to the free-fall velocity Constant temperature and density are 
set at the outer boundary to denote the nebula conditions. In order to approximate the 
effects of the tidal potential, we added a simple prescribed outward gravitational force in 
the equation of motion (Equation [19]). We can then adjust the location of the Roche radius 
where gravity vanishes in order to assess the effects of this tidal barrier. 

For comparison, we first constructed a test model for standard Bondi accretion with 
the nebula conditions of T ne b ~ 50K and p ne b ~ 1 x 10~ 14 g/cm 3 (Model in Table [p. 
These parameters are relevant in the outer (> 5 — 10 AU) regions after a substantial fraction 
(~ 10~ 3 — 10~ 2 ) of the minimum mass nebula has already been depleted through global 
depletion or gap formation. In this, we neglect the radiative energy transfer and tidal barrier 
effects so that we can compare our numerical results with the the analytic Bondi solution. 
We find a gas accretion rate of around 3 x 10 17 g/s from our numerical calculations, which 
i s in agreement wi th the analytic Bondi solution M ~ 1.4 x 10 n ( ^) 2 (j^r)( 1 Q a k < ^ D J-i )^ 3 g/s 



flFrank et alJll992h . 



We constructed another three models to investigate the tidal barrier effects, with differ- 
ent parameters listed in Tabled] (Model 1, 2, 3). In the table, T neb denotes the temperature of 
the nebula gas, which corresponds to the sound-speed at the outer boundary of the envelope. 
The third column shows the ratio of Roche radius to the Bondi radius. The last column lists 
the average gas accretion rate for each model. The outer boundary condition for density is 
1 x 10~ u g/cm 3 for all the models. 
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Model No. 


T ncb (K) 


Rr/Rb 


M (g/s) 





50 




3.0 x 10 17 


1 


50 


0.5 


1.6 x 10 14 


2 


50 


1.0 


1.7 x 10 17 


3 


100 


1.0 


6.2 x 10 16 



Table 1: Parameters for the single, one-dimensional models. Model shows Bondi accretion, 
which does not include radiative energy transfer or the tidal barrier. T ne b denotes the tem- 
perature of the nebula gas, while M is the average gas accretion rate calculated from the 
model. 
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We evolved all the models to steady-state, where the accretion rate became constant 
with radius. Figure (fT4"j) shows the temperature profiles of the three models after they have 
reached steady-state. As we can see from the figure, the outer regions of the envelopes are 
optically thin, cool, and isothermal for all the models. Figure (fT5"j) plots the density profiles 
of the three models. The model with Rr/R b ~ 0.5 shows a much larger density gradient and 
thus much smaller density at the location of Roche radius than the other two models. At 
small radii, the pressure gradient does not significantly perturb the nearly free-falling motion 
of the gas. However, as the stellar gravity acts to repel the gas outside the Roche radius, 
sufficiently large positive pressure and density gradients are needed to overcome the tidal 
barrier. Therefore, in the case of Rb > Rr, this gradient leads to a density at R R that is much 
smaller than that in the disk. Since free fall onto the core can only proceed inside Rr and 
the flow across Rr is limited by the speed of sound, the mass accretion rate is much reduced 
from the conventional Bondi formula. Figure ffToT) shows the evolution of gas accretion rate 
with time at the locations of Roche radius for three models. Nearly constant accretion rate 
indicates that all the models have reached steady state after approximately 6 x 10 9 s. Since 
the flow between R B and Rr is significantly perturbed by the gas pressure gradient and 
softened gravity, this time scale must be at least several times the local dynamical free-fall 
time scale (~ 5 — 8 x 10 8 s). However, in contrast to the 2D simulations, the flow does not 
need to adjust over the liberation period of flows in the horseshoe region and the synodic 
period of the flow near the separatrix at the Lagrangian points. 

From the two models with the same Rr/ Rr but different nebula temperature (Model 2 
and Model 3), it's found that higher nebula temperature leads to slightly lower gas accretion 
rate. This behavior is similar to what we find in the Bondi solution, in which Moc c7 3 (oo). 
When comparing the models with the same nebula temperature but different Rr/ Rr (Model 
1 and Model 2), we find that the accretion rate is reduced by three orders of magnitudes for 
the model with R B > R R . Thus, the tidal barrier significantly suppresses in the gas accretion 
across the Roche radius. Even if the models are surrounded by nebula gas with the same 
density, only a small fraction of that gas can be accreted by the model with R B > Rr. 
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Fig. 14. — Temperature profiles as a function of radius for ID models 1, 2, and 3 after they 
have reached steady-state. The parameters for each model are listed in Table (CQ). 
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Fig. 15. — Density profiles as a function of radius for the three ID models after they have 
reached steady-state. The denotations are the same as in Figure (j!4p . 
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Fig. 16. — Gas accretion rate as function of time for the three models. Different line types 
correspond to the different models listed in Table (Q. 
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We further compute two series of disk models (4 and 5) for planetary masses ranging from 
30M e to 100 and 300M e , respectively. Models in series 4 have (T, p neb ) of (100K, 10" 12 g cm" 3 ), 
while models in series 5 have (T,p neb ) of (50K, 10~ 12 g cm -3 ). We determine the magnitude 
of M, as a function of M p after it has attained a steady state value at 10 AU (see Figure [T7|) . 
These results are in good agreement with the relationship in Equation (|23|) . Based on these 
results, we determine the total mass of the protoplanet after 20 Myr of accretion onto a core 
with an initial mass of 30M ffi (see Figure ITS]) . In both model series we neglect the decline 
of gas density associated with formation of the gap in the disk, but include that due to the 
global depletion over a time scale Td ep = 3Myr. Despite these simplifications, the asymptotic 
values of M p (100 and 300 M e for T ne & = 50K and 100K respectively) are comparable to 
that found in Equation fl33l with the appropriate boundary conditions. 



4. Summary and Discussion 

In this paper, we present evidence to suggest that the asymptotic mass of protoplanets 
is determined by the structure of their nascent disk. The main physical process involved here 
is that the tidal potential of the central star provides a barrier to the gas flow in the vicinity 
of the protoplanet. Through its tidal interaction with the disk, a growing protoplanet exerts 
an increasingly strong torque which first leads to the formation of a gap. Although the disk 
gas may continue to diffuse into the gap, it is forced to flow along the horseshoe stream 
lines. In the absence of any dissipation, the vortensity of the flow along the stream lines is 
conserved. But near the protoplanets, shock dissipation can lead to vortensity generation 
and dissipation. 

In the limit that the protoplanets' mass is small and its Roche radius is smaller than 
the disk thickness, the tidal potential barrier is shallow and modest surface density variation 
would provide adequate pressure gradient to overcome the tidal potential barrier while pre- 
serving vortensity along the stream line. However, in cold disks with thicknesses less than 
the protoplanets' Roche radius, a very large surface density gradient is required to create 
the critical pressure gradient force needed to overcome the tidal barrier along the stream 
line. Consequently, mass supply from the feeding zone into the vicinity of the protoplanet 
is quenched. In a pressure-free gas disk, the horseshoe stream lines of cold gas cannot join 
onto those in the disk around the protoplanet. Only particles with a limited range of initial 
orbital parameters that do not follow streamlines can directly strike the protoplanet. 



Based on this consideration, we suggest that the asymptotic mass of the protoplanet 
is determined by t he condition Rr ~ H. T his condition is similar to the thermal criterion 
for gap formation flLin Papaloizoul Il986l ). When the planet approaches this value, the 
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Fig. 17. — The steady-state gas accretion rate as a function of M p for models series 4 and 
5 are represented by solid and dotted lines respectively. Both simulations show a significant 
decrease in M, but it occurs at much lower mass for the low temperature nebula. 
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Fig. 18. — The planet's mass after 20 Myr of gas accretion, starting with M p = 30M e 
initially. The gas depletion time scale is set to be r dep = 3 Myr. The solid line designates 
the model series with T ne {, = 50K, while the dotted line is for the model series with T ne & = 
lOOf^. As suggested from Figure [T71 and seen in Figured the asymptotic mass of the lower 
temperature model series is significantly lower then in the high T ne b model series. 
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mass supply into th e gap is much reduced and the dynamics of the gas is slightly altered 
( IBryden et al.lll999l ). Although these consequence do not directly terminate protoplanets' 
growth, they contribute to the reduction of the accretion rate onto the protoplanet. 

The gaseous material which enters into the protoplanet 's Roche lobe forms a circum- 
planetary disk. Angular momentum is transferred from the circumplanetary disk by both 
tidally induced spiral shocks and turbulent viscosity. Initially, the flux through the circum- 
planetary disk is hi gh and the disk environment is prohibitive for the formation of regu- 
lar satellite systems (ILunine fc Stevenson! Il982l ; ICanup fc Ward! 120021 ; iMosqueira fc Estrada 



20031 : ICanup fc Ward! l2006f ). However, our results indicate that during the final stage of 
the protoplanet 's growth, the characteristic timescale for protoplanetary accretion increases 
beyond the circumstellar disk depletion timescale. It is at this stage that p hysical environ- 
ment in the circumplanetary disks become less hostile to satellite formation (ICanup fc Ward 
20061 ). The studies of satellite formation process is beyond the scope of the present paper 
and they will be presented elsewhere. 

In the final analysis, we are interested in the origin of the present-day mass distribution 
of extrasolar planets whic h, is observe d to be dN/dM cx M~ l to M~ 1A in the range of 



a fraction to several Mj (IMarcy et al.l 120051 ). This distribution also has a cut-off below 



M p ~ 10 — 20Mj. Since the host stars of most known extrasolar planets have masses 

1 to q~ 1A in 



M* ~ M Q , the corresponding q 
the range ~ 10~ 3 . 



Mp/M* has a similar distribution dN/dq ~ q 



Scaling the formula derived in §2.31 to the values used in this simulation and assuming 
a dissipation rate of T^ ep = 3 — 5Myr and /3 = 1, M p is comparable to the mass of Saturn for 
an aspect ratio of 0.04 and that of Jupiter for an aspect ratio of 0.07. T he most favorabl e 
location for the formation of gas giant planets is near the snow line a; ce (jlda fc Linll2004l ). 
The magnitude of aj ce is determined by both the gas accretion rate within the disk and the 
intensity of stellar irradiation and dur ing the classical T-T auri phase, and evolve over 

nearly an order of magnitude in scale (IGaraud &: Linll2006l ). The disk aspect ratio at the snow 
line (H (a ice ) /a ice ) can also change by a factor of two. Depending on the epoch of gas-giant 
planet formation, the variation in (H (a ice ) /a ice ) can lead to an order of magnitude spread 
in the asymptotic mass of the planets. Thus, the observed planetary mass distribution may 
be used as to infer the ratio of the planet formation and disk depletion timescales. Detailed 
simulations of such models will be presented elsewhere. 

Target lists for planet searches using the radial velocity technique have recently been 
expanded to include M-dwarf stars. However, the detection frequency for Jupiter mass 
planets around these low-mass stars appears to be much lower than that around G dwarfs. 
Until recently, the only known member of this population with known Jupiter-mass planets 
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was GJ 876. The lack of massive gas-giants is also conspicuous in the microlensing survey 
for extrasolar planets (IBeaulieu et al.ll2006l ). The difficulty for gas-giant planet formation 
may be due to the in adequate amount of heavy elements for assembling sufficiently massive 
protoplanetary cores (llda fc Linl 120041 ) or inability to overcome type I migration and r etain 
sufficiently mass cores to initiate the process of rapid gas accretion (jlda fc Linl 120061 ) . In 
this paper, we show that even in the low-probability event that proto gas-giant planets 
may form, their asymptotic mass may severely limited by the tidal barrier in accordance 
with Equation (1551) . If the surface densities of both gas and dust increases with Ml as 
suggested by observationally inferred M's from the disks onto their central stars, the onset 
of dynamical gas accretion would occur at much later epoch with a relatively small total gas 
content around M dwarfs than around solar type stars. At this stage, the value of H(a lce )/a 
is also an slowly increasing function of M*. All of these effects combine sensitively to quench 
the gas accretion rate and severely limit the asymptotic mass of the planets well below that 
of Jupiter. 

The co-existence of multiple gas-giant planets around the same host star suggests that 
their asymptotic mass is limited by local tidal truncation (through gap formation) rather 
than the global depletion of the disk. The asymptotic mass distribution within such systems 
is determined by both the planets' formation location and epoch. In a steady-state disk, 
the value of the aspect ratio H/a de termined from the disk mid-plane structure generally 
increases with a (IGaraud fc Linl 120061 ) . In a viscous evolving disk, E in the outer regions of 
the disk decreases with a much more rapidly than in a steady disk. In general, H/a attains 
a maximum at a radial location a max which is an increasing function of the the accretion 
rate through the disk. During the depletion of the disk, a max decreases. Outside a max , the 
disk is shielded from the stellar irradiation and the density scale decline rapidly with radius. 
This process may determine the mass distribution within multiple gas giant planet systems. 
It may also limit the domain of gas-giant formation. 
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